
subroutine InitEphem(Ephem, PathToEphem)
    use calc_eph
    implicit none
    integer Ephem
    character*80 PathToEphem
    character*80 filename
    
    character*80 EXTNAMFIL/""/
    COMMON/ETHEM_FILE_NAME/EXTNAMFIL
    
    logical:: EpmUsedAsEphem
    COMMON /EPMorDEused/EpmUsedAsEphem
    !--------------------------------
    
    EpmUsedAsEphem = .FALSE.
    if(Ephem .EQ. 1) then 
        filename = 'epm2008\'
    else if(Ephem .EQ. 2) then ! de403
        filename = 'de403.dos'
    else if(Ephem .EQ. 3) then ! de405
        filename = 'de405.dos'
    else if(Ephem .EQ. 4) then ! de414
        filename = 'de414.dos'
    else if(Ephem .EQ. 5) then ! de421
        filename = 'de421.dos'        
    else if(Ephem .EQ. 6) then ! de423
        filename = 'de423.dos'        
    else
        write(*,*) 'Unknown Ephem number', Ephem
        return 
    end if
    
    call strcat(EXTNAMFIL, PathToEphem, filename)
    write(*,*) EXTNAMFIL
    
    if(Ephem .EQ. 1) then
        EpmUsedAsEphem = .TRUE.
        call InitEpm(.FALSE., EXTNAMFIL)
        CALL initGMsEPM()
        CALL InitEphemTimeRangeEPM()
    else !de
        CALL initGMsDE()
        CALL InitEphemTimeRangeDE()
    end if
end

! =========================================================

subroutine strcat(res,str1,str2)
       implicit none
       character*(*) res,str1,str2
       integer*4 len1,len2
       res=trim('')
       len1=len_trim(str1)
       len2=len_trim(str2)
       res(1:len1) = str1(1:len1)
       res(len1+1:len1+len2+1) = str2(1:len2)
       res(len1+len2+3:len1+len2+3) = " "
end

! =========================================================

subroutine GetCoordVelRegBody(TM, RegBodyX, RegBodyV, CenterBody)
    use calc_eph
    implicit none
    logical:: status         ! epm success
    REAL*8 EthenTimeRange(3)
    COMMON /EthenTimeRange/EthenTimeRange
    logical:: EpmUsedAsEphem
    COMMON /EPMorDEused/EpmUsedAsEphem
    
    CHARACTER*10 name_body(11)
    COMMON /name_body/ name_body
    
    REAL*8 RegBodyX(3,11), RegBodyV(3,11), epmXX(3), epmVV(3)
    REAL*8 TM
    REAL*8 R(6)
    INTEGER i, j
    integer CenterBody
    integer InnerCenterBoby
    !---------------------------------    
    
    if( (TM .LT. EthenTimeRange(1)) .OR. (TM .GT. EthenTimeRange(2)) ) then
        WRITE(*,*) '~~~~',TM, ' out of EPHEM time range ',&
            EthenTimeRange(1), EthenTimeRange(2)
        return
    end if
    
    InnerCenterBoby = CenterBody
    if(InnerCenterBoby .EQ. 0)    InnerCenterBoby = 12
        
    do i=1, 11
        if(EpmUsedAsEphem) then
            call calc_EPM(i, InnerCenterBoby, 0, TM, RegBodyX(1, i), RegBodyV(1, i), status)
            if (.NOT.status) then
                write(*, *) i, ' date error '
                return
            end if  
        else
            CALL  Pleph ( TM, i, InnerCenterBoby, R )
            IF (1.eq.0) then
                write(*,*) name_body(i)
                write(*,*) 'pos: ', (R(j),j=1,3)
                write(*,*) 'vel: ', (R(j),j=4,6)
                write(*,*) ''
            end if
            do j=1,3
                RegBodyX(j, i) = R(j)
                RegBodyV(j, i) = R(j+3)
            end do
        end if
    end do
end

! =========================================================

subroutine GetCoordOnlyRegBody1(TM, RegBodyX, CenterBody)
    implicit none
    REAL*8 TM
    REAL*8 RegBodyX(3,11), RegBodyV(3,11)
    integer CenterBody
    call GetCoordVelRegBody(TM, RegBodyX, RegBodyV, CenterBody)
end

! =========================================================

subroutine GetCoordOnlyRegBody(TM, RegBodyX, CenterBody)
    use calc_eph
    implicit none
    
    REAL*8 EthenTimeRange(3)
    COMMON /EthenTimeRange/EthenTimeRange
    logical:: EpmUsedAsEphem
    COMMON /EPMorDEused/EpmUsedAsEphem
    
    integer CenterBody
    integer InnerCenterBoby
    logical:: status         ! epm success
    
    REAL*8 RegBodyX(3,11), RegBodyV(3)
    REAL*8 TM
    REAL*8 R(6)
    INTEGER i, j
    
    INTEGER LIST(12),IPT(3,13), NUMDE, NCON
    LOGICAL KM, BARY/.TRUE./
    REAL*8 PVSUN(6), ET2(2), RRD(6), PV(6, 12), SS(3),CVAL(400), AU, EMRAT
    COMMON/STCOMX/KM,BARY,PVSUN
    COMMON/EPHHDR/CVAL,SS,AU,EMRAT,NUMDE,NCON,IPT
    
    !---------------------------------    
    if( (TM .LT. EthenTimeRange(1)) .OR. (TM .GT. EthenTimeRange(2)) ) then
        WRITE(*,*) '~~~~',TM, ' out of ETHEM time range ',&
            EthenTimeRange(1), EthenTimeRange(2)
        return
    end if
    
    InnerCenterBoby = CenterBody
    if(InnerCenterBoby .EQ. 0)    InnerCenterBoby = 12
    
    if(EpmUsedAsEphem) then
        do i=1, 11
            call calc_EPM(i, InnerCenterBoby, 0, TM, RegBodyX(1, i), RegBodyV, status)
            if (.NOT.status) then
                write(*, *) i, ' date error '
                return
            end if  
        end do
    else 
        DO I=1,10
            LIST(I)=1
        ENDDO
        DO I=11,12
            LIST(I)=0
        ENDDO
        
        ET2(1) = TM
        ET2(2) = 0.D0
        
        if(InnerCenterBoby .NE. 12) BARY=.FALSE.
        CALL STATE(ET2,LIST,PV,RRD)
        
        do i=1, 9
            do j=1,3
                RegBodyX(j,i)=PV(j,i)
            end do
        end do
        do i=1,3
            RegBodyX(i,3) = PV(i,3) - PV(i,10)/(EMRAT+1.0d0)  !EARTH
            RegBodyX(i,10) = RegBodyX(i,3) + PV(i,10)  !MOON
        end do
        do i=1,3
            RegBodyX(i,11)=PVSUN(i)  ! SUN
        end do

        endif
end

! =========================================================

subroutine FixBodyPos(CenterBody, TM, X, V, toBary)
    implicit none    
    
    LOGICAL toBary
    REAL*8 TM
    REAL*8 X(3), V(3)
    REAL*8 RegBodyX(3,11), RegBodyV(3,11)
    INTEGER CenterBody, j

    !--------------------------------
        
    if ( (CenterBody .NE. 0) .AND. (CenterBody .NE. 12) ) then
        call GetCoordVelRegBody(TM, RegBodyX, RegBodyV, 0)
        if (toBary) then
            write(*,*) 'place body into barycenter'
            do j = 1, 3
                X(j) = X(j) + RegBodyX(j, CenterBody)
                V(j) = V(j) + RegBodyV(j, CenterBody) 
            end do
        else
            write(*,*) 'place body into heliocentric'
            do j = 1, 3
                X(j) = X(j) - RegBodyX(j, CenterBody)
                V(j) = V(j) - RegBodyV(j, CenterBody) 
            end do
        end if
                
        write(*,*) 'Sun pos:'
        write(*,*) (RegBodyX(j, CenterBody), j=1,3)
        write(*,*) (RegBodyV(j, CenterBody), j=1,3)
        write(*,*) 'New Body pos:'
        write(*,*) (X(j), j=1,3)
        write(*,*) (V(j), j=1,3)
    end if
end
